      program id5_L0
      implicit real*8 (a-h,o-z)
      parameter(ndec=12,npat=4096,ngrid=51)
c
      dimension pa(ndec,2),tl(npat+1,2),pz(ngrid),g(ngrid),bgrid(ngrid)
      dimension er1(ngrid,ngrid),er2(ngrid,ngrid)
      dimension in(npat,ndec)
      common /pat/in
c
      write(*,'('' ***** id5_L0 *****'',
     &//"    Level-0 versus AA"/)')
c
      open(unit=10,file='pat4096.data')
      do 10 n=1,npat
         read(10,*)idum,(in(n,k),k=1,ndec)
c         write(*,'(i5,x,12i2)')n,(in(n,k),k=1,ndec)
   10 continue
c
      do 15 i=1,ngrid
         bgrid(i) = 0.02*(i-1)
         k = ngrid + 1 - i
         if (k.eq.ngrid) then                                
            pz(i) = 0.999d0
         else
            pz(i) = 0.5d0 + 0.01d0*(k-1)
         endif
         g(i) = 2.d0*dlog(pz(i)/(1.d0-pz(i)))
   15 continue
c
      x = dlog(0.5d0)*float(ndec)     !LL of L0 model
      do 16 j=1,ngrid
            er1(ngrid,j) = 0.5d0                    ! gam = 0
            er2(ngrid,j) = 0.5d0                    
   16 continue
      do 100 i=1,ngrid-1
         gam = g(i)
         do 50 j=1,ngrid
            beta = bgrid(j)
c
            call fpx(gam,beta,pa)
            call fl(pa,tl)
            call sort(npat+1,tl)
c
            p0 = 1.d0/float(npat)
            f0 = 1.d0
            fa = 0.d0
            do 30 n0 = 1,npat+1
               if (tl(n0,1) .ge. x) then
c                 write(*,'(3g12.4)')x,f0,fa
                 goto 31
               endif
               f0 = f0 - p0
               fa = fa + tl(n0,2)
   30       continue
   31       er1(i,j) = fa                    ! Prob(LL<x|Ha)
            er2(i,j) = f0                    ! Prob(LL>x|H0)
   50     continue
          if (er1(i,j).gt.er2(i,j)) er2(i,j)=er1(i,j)
  100 continue
c
      write(*,'(//"Max of Type I & II Errors:"/)')
      write(*,'(8x,51f8.4)')(bgrid(j),j=1,ngrid)
      do 110 i=1,ngrid
         write(*,'(52f8.4)')pz(i),(er2(i,j),j=1,ngrid)
  110 continue
c
      end
c***********************************************************************
      subroutine fpx(gam,beta,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12)
      dimension px(ndec,2)
c
      y1 = 10.d0
      y2 = 12.d0
      y3 = 15.d0
      y4 = 5.d0
      y5 = 6.d0
      y6 = 8.d0
      pjb = (1.d0+beta)/3.d0
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(1,1) = za/(za+zb)
      px(1,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y1
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(2,1) = za/(za+zb)
      px(2,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(3,1) = za/(za+zb)
      px(3,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y2
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(4,1) = za/(za+zb)
      px(4,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(5,1) = za/(za+zb)
      px(5,2) = zb/(za+zb)
c
      za = 0.5d0*y1
      zb = 0.5d0*beta*y3
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(6,1) = za/(za+zb)
      px(6,2) = zb/(za+zb)
c
      za = y1*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(7,1) = za/(za+zb)
      px(7,2) = zb/(za+zb)
c
      za = y2*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(8,1) = za/(za+zb)
      px(8,2) = zb/(za+zb)
c
      za = y3*beta*1.d0/3.d0
      zb = y1*1.d0/3.d0
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(9,1) = za/(za+zb)
      px(9,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y4*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(10,1) = za/(za+zb)
      px(10,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y5*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(11,1) = za/(za+zb)
      px(11,2) = zb/(za+zb)
c
      za = y4*2.d0/3.d0
      zb = y6*pjb
      za = dexp(gam*za)
      zb = dexp(gam*zb)
      px(12,1) = za/(za+zb)
      px(12,2) = zb/(za+zb)

c      do 10 k=1,ndec
c         write(*,'(i2,2g12.4)')k,(px(k,j),j=1,2)
c   10 continue
c
      return
      end
c***********************************************************************
      subroutine fl(px,tl)
      implicit real*8 (a-h,o-z)
      parameter (ndec=12,npat=4096)
      dimension px(ndec,2),tl(npat+1,2)
      dimension in(npat,ndec)
      common /pat/in
c
      pt = 0.d0
      do 50 n=1,npat
         t1 = 0.d0
         do 10 k=1,ndec
            if (in(n,k) .eq. 1) then
               t1 = t1 + dlog(px(k,1))
            else
               t1 = t1 + dlog(px(k,2))
            endif
   10    continue
         tl(n,1) = t1
         tl(n,2) = dexp(t1)
         pt = pt + tl(n,2)
c      write(*,'(i5,2x,3g12.4)')n,t1,tl(n,2),pt
   50 continue
      tl(npat+1,1) = 1.d2
      tl(npat+1,2) = 0.d0
c      write(*,'("pt = ",g12.5)')pt
c      stop
c
      return
      end
c***********************************************************************
      subroutine sort(n,tl)
      implicit real*8 (a-h,o-z)
      dimension tl(n,2),a(n),b(n),indx(n)
c
      do 5 i=1,n
         a(i) = tl(i,1)
         b(i) = tl(i,2)
    5 continue
      call indexx(n,a,indx)
c
      do 10 i=1,n
         j = indx(i)
         tl(i,1) = a(j)
         tl(i,2) = b(j)
c      write(*,'(i5,x,2g12.4)')i,tl(i,1),tl(i,2)
   10 continue
      return
      end
c***********************************************************************
      SUBROUTINE INDEXX(N,ARRIN,INDX)
      implicit real*8 (a-h,o-z)
      DIMENSION ARRIN(N),INDX(N)
      DO 11 J=1,N
        INDX(J)=J
11    CONTINUE
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          INDXT=INDX(L)
          Q=ARRIN(INDXT)
        ELSE
          INDXT=INDX(IR)
          Q=ARRIN(INDXT)
          INDX(IR)=INDX(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            INDX(1)=INDXT
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
          ENDIF
          IF(Q.LT.ARRIN(INDX(J)))THEN
            INDX(I)=INDX(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        INDX(I)=INDXT
      GO TO 10
      END


